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The open hole tensile and compressive strengths are important design 
parameters in qualifying fiber reinforced laminates for a wide variety of 
structural applications in the aerospace industry. In this paper, we present 
a unified model that can be used for predicting both these strengths (tensile 
and compressive) using the same set of coupon level, material property data. 
As a prelude to the unified computational model that follows, simplified ap- 
proaches, referred to as “zeroth order”, “first order”, etc. with increasing 
levels of fidelity are first presented. The results and methods presented are 
practical and validated against experimental data. They serve as an intro- 
ductory step in establishing a virtual building block, bottom-up approach 
to designing future airframe structures with composite materials. The re- 
sults are useful for aerospace design engineers, particularly those that deal 
with airframe design. 


I. Introduction 

Predicting the open hole compression strength (OHCS) and the open hole tensile strength 
(OHTS) of a laminated fiber reinforced composite is an important consideration in material 
qualification for aerospace structural design. Experimental results for fiber reinforced com- 
posite panels with holes and cut-outs show a significant drop in compressive strength as a 
function of hole size due to stress gradients generated by the hole and because of the large 
stresses that concentrate at the hole edge. The situation is similar for tensile loading. 

In traditional multi-layered laminates, the failure in compression is generally dominated 
by the 0° ply kink banding mechanism 1, while in tension, transverse cracks in the 90 
degree layers, and matrix microcracks between fibers in off-axis layers are observed prior 
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to the onset of catastrophic failure, the latter resulting due to failure of the significant 
load-bearing zero plies, and sometimes also accompanied by delamination. Even though 
the problem of predicting compressive strength of laminated composite panels with a cut- 
out is complex, it can be simplified by focussing attention on the failure of the 0° lamina 
within the laminate. There are two main points that need to be taken into consideration; 
first the compressive strength of the 0° lamina and second, the influence of the hole on the 
stress distribution around the hole. Prior studies, 1 ’ 2,4,5 have examined the problem in detail 
and have developed both analytical and numerical methods based on careful experimental 
results to predict OHCS with varying degrees of sophistication in the modeling. In tension 
loading, in addition to transverse cracking in the 90 degree layers, matrix microcracking and 
fiber-matrix splitting in off-axis layers, and delaminations can all influence tensile strength. 
Hence, both damage and failure have to be taken into consideration in modeling for OHTS 
prediction. 

At the zeroth order, we present a simple engineering approach using a combination of the 
Considere construction, 7 and Lckhnitskii’s 8 formulation to approximate OHCS. For OHTS, 
we use a modified Lckhnitskii’s analysis which uses a strain based failure index for predicting 
strength. At the first order, these models are refined by a combination of micromechanics 
and classical lamination theory, implemented using the finite element method to capture 
kink-band formation using a global-local approach, 5 for accurate prediction of OHCS for a 
class of popular laminates. For OHTS, a similar modeling strategy with the smeared crack 
approach (SCA), described in Heinrich and Waas, 13 can be used for a single fiber reinforced 
lamina to predict OHTS. For both cases, the modeling uses constitutive properties of the 
fiber and in-situ lamina shear-stress strain response. In the second order analysis, which has 
higher fidelity, we employ the thermodynamically-based Schapery work potential damage 
theory for capturing the in-plane damage in each lamina, and we also use the enhanced 
Schapery theory (EST), as described in Pineda et al, 15 to model in-plane failure including 
post-peak softening using the crack band model. 22 In this paper, results for OHCS are 
predicted using the zeroth, first and second order analysis methods, while results for OHTS 
are predicted using the zeroth order and second order analysis methods. 

II. Zeroth order analysis 

In multi-directional laminates, compression strength is generally dominated by the 0° 
ply kink banding failure. Even though the problem of predicting compressive strength of 
laminated composite with a cut-out is complex, it can be simplified by focusing attention 
on the failure of 0° laminae within the laminate. More specifically, the failure of 0° laminae 
occurs at the location of maximum stress, which, in case of laminates with cut-outs like 
a hole, is at the hole edge. Hence, the simplest approach would involve; (a) knowing the 
compressive strength of 0° lamina and, (b) knowing the value of stress at the hole edge. 
Therefore, the compressive strength of laminates with a hole is taken to be the applied stress 
which causes the maximum stress at the hole edge to reach the compressive strength of the 0° 
lamina. This approach is implemented using a combination of the Considere construction, 7 
and Lekhnitskii’s 8 formulation. Considere construction is a geometric construction which 
utilizes the lamina shear stress-strain response to obtain the compressive strength of a 0° 
lamina (X?) for a given fiber misalignment q i, as shown in figure 2. 
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where, (f> is the fiber misalignment, 772 is the shear stress and 712 shear strain. The compres- 
sive strength of a laminate, X L C , is obtained by, 


Vf = V»(||) 


( 2 ) 


where, E\\ and E\\ are the elastic moduli of the laminate and the 0° ply, respectively, in the 
direction of compressive load. Further, is the compressive strength of a 0° lamina for a 
given fiber misalignment <fi. 

The equation above provides the compressive strength of a laminate without any stress- 
raisers (such as at cutout edges). In the presence of a hole, the maximum stresses occur at 
the hole edge due to non-uniform stress distribution. Using the Lekhnitskii 8 solution, Filiou 
and Soutis,, 11 have provided an approximation to the stress distribution in an orthotropic 
infinite plate with a hole (shown in figure (3)). The equation below provides an accurate 
description of the axial stress, (a yy (x, 0)) distribution due to a far held stress, cr^, 
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Using the above equation, we can approximate the stress at the edge of the hole (for 
x = R) as, 


CyyjR, 0 ) 

&OO 


~3 H a 


( 4 ) 


Assuming the material at the edge of the hole, where the stress is maximum, fails at the 
strength provided by equation 4, above, ie; <J yy (R , 0) ~ X T C J we can obtain the corresponding 
far held stress. This far held stress will be the critical strength of the laminate with a hole 
in compression, which corresponds to the zeroth order approximation. 


L ^ a yy{ 0, R) _ X® . Ejj 

c ~ 3 H a 3 HaE q J ( ’ 

Thus, this equation provides the OHCS of a laminate (X C L , in terms of the pristine lamina 
compressive strength, (A( J ). The approach presented here can also be used to estimate the 
open hole compressive strength for other cutout shapes and for multiaxial planar loading, 
provided an expression that relates the maximum stress at the cutout edge in terms of the 
applied far-held loading is available. 

To estimate OHTS, we utilize a similar approach except we assume that the laminate 
will fail when the strain in the 0° ply reaches its ultimate tensile strain (X® t ) obtained from 
a coupon [O4] tensile test. A uniform through-the-thickness strain is also assumed in this 
model. Equation 5 will now be modified to, 
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Equations 5-6 provide the zeroth order estimates for OHCS and OHTS. With the as- 
sumptions used for this simplified analysis, there is a significant drawback. The model does 
not account for material degradation due to in-plane damage and inter-laminar failure. In- 
plane damage and delamination is seen to be negligible in laminates with high percentage 
of 0° plies, but that is not always the case with laminates with higher percentage of off-axis 
plies, such as quasi-isotropic laminates. This method should only be used as a starting point 
for the more refined methods that will be next described. 


III. First order analysis 

First order analysis is based on the same mechanism as the zeroth order approach, that 
the failure in laminates under compressive loading is due to kink banding. However, instead 
of analyzing the kink banding of a pristine 0° lamina, and the hetrogeneous stress held 
due to a cut-out separately, the kink band formation in a 0° lamina in a hetrogeneous 
stress held is analyzed. This is achieved through a micro-mechanical, finite element model 
of the 0° lamina with a cut-out. In the micro-mechanics model, the hbers are explicitly 
modeled as homogeneous orthotropic continua and the matrix is modeled as an isotropic 
material with in-situ elastic-plastic properties. Fibers are misaligned by 1° to account for 
average hber misalignment angles commonly observed in carbon composite materials. 3 This 
misalignment is the same as the one used in the zeroth level analysis. The micro-mechanics 
model, however, is a computationally expensive method, because of the difference in length 
scale of hber and lamina. To model a complete lamina using micro-mechanics will require 
thousands of hber-matrix layers, which is computationally prohibitive. To overcome the 
issue of computational cost the model is reduced by analysing a micro-region around the 
hole. Boundary conditions to be applied on the micro-region are obtained from a global 
homogeneous model. This method is called the global-local approach as defined in the work 
by Ahn and Waas 5 and Davidson and Waas. 6 In some reports, this method has also been 
referred to as the embedded cell method (see, for example, Llorca et al 19 ). 

The local-global approach implemented here is a three step procedure (each subsequent 
step relies on the results of the prior step) using a macro-mechanical model and this is as 
depicted in figure (4). Each step is an individual simulation and only the boundary infor- 
mation is passed from one step to the other. The macro-model, outside the micromechanics 
sub-region, is a 2D half symmetric model with homogenized transversely isotropic lamina 
properties of the 0° lamina. Step 1, is a linear elastic step where a known far-held displace- 
ment is applied in the axial (y) direction and the displacements (A X) y) at nodes corresponding 
to the micro-region boundary (A — B — C — D) are obtained. These displacements A XjV are 
then applied on the isolated local micro-mechanics model boundary and another linear elas- 
tic analysis is performed in step 2, henceforth called ’’Local Linear Model”. In step 3, the 
reaction forces, (R x , y ), obtained in step 2 are applied on boundaries ( A — D)&,(B — C) along 
with A XjV on (A — B)!k(C — D) and a non-linear Riks response analysis is conducted (’’Local 
Non-linear Model”) on the isolated local micro-mechanics model. Only step 3 is an incre- 
mental analysis while steps 1 and 2 are linear and hence done only once. Step 2 is needed 
to obtain the consistent reaction forces on parts of the boundary of the micro-region, since 
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step 3 uses mixed boundary conditions on the micro-region boundary. 

Typically, as seen in past studies, there are no boundary conditions imposed on the free 
edge of a local model. This however, can lead to incorrect results because the local model, 
which is hetrogeneous, will deform in a manner different from the global model. In figure 5, 
the un-deformed and deformed profiles of global and local model are shown. Pinching of the 
corners of the micro-region is observed as compared to the global model, during deformation, 
which will lead to artificial stress concentration developing at the corners. These stress 
concentrations will cause premature kink formation at the corners. Hence, application of 
reaction forces at the free boundaries is important to ensure that the local model follows 
the global deformation profile, in the elastic regime. By applying reaction forces instead of 
displacements, the left and right boundaries are kept free of constraints, thereby allowing 
for the kink band to form without displacements being constrained on the boundaries. 

The axial reaction force R y obtained from steps 1 & 2 provide the linear relation between 
the local and global model as a simple ratio R = (R9 lobal / R l ° cal ) . Now, from step 3, we 
get the micro-kink band initiation load (R l k cal ) which corresponds to the global compressive 
strength of the zero lamina due to kink banding, X jy as ; 


X° N = (R 9 k lobal ) I A = (R R[ ocal ) /A (7) 

In the above expression, A, is the load bearing cross-sectional area of the specimen in 
the far held. The size of the micro (local) region is an unknown of the global-local method. 
A converged solution is obtained when the peak compressive strength of the first order 
analysis shows no difference with respect to the size of the micro-region. The size of the 
region was scaled based on the ratio of width (length A-B) to the height (length B-C) of the 
micro-region. 

Typical response of the micro-region is shown in figure (6) and corresponding snap- 
shots at different stages of loading leading to kink banding failure are provided. Figure (8) 
shows the scaling study to determine the converged micro-region size and the corresponding 
compressive strength of the 0° ply. The results are normalized by the compressive strength 
of the 0° ply with 0 = 1° obtained from the Considere construction ie; the un-notched zero 
ply strength. Using the converged value of compressive strength from the 2D global local 
analysis, the laminate compressive strength can now be obtained by simple scaling, 


rpL 

x^O / ^11 

AV 


( 8 ) 


As with zeroth order, this first order approach provides compressive strength prediction 
for laminates that show kink band formation in the 0° lamina. Hence, this method of 
predicting compressive strength is predicated on zero ply kink banding being the dominant 
failure mechanism. However, this may not provide satisfactory predictions in cases where 
off-axis laminae dominate, ie; in cases where the percentage of 0° laminae is small and where 
there is significant non-linear behaviour before failure, as in pure angle-ply laminates or 
laminates with a large percentage of angle plies. 
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IV. Second order analysis 


Zeroth and first order analysis methods do not account for stiffness degradation of off-axis 
laminae. The difficulty in analytically modeling off-axis plies is because significant damage 
accumulates in these plies during loading. A second order modeling method should capture 
both damage and failure. This is achieved by modeling damage using Schapery theory and 
failure using the crack-band model ,, 22 which is referred to as Enhanced Schapery Theory 
(EST). 

The EST formulation developed by Pineda et al 15 is used, which assumes that three major 
intralaminar mechanisms are responsible for all observed non-linearities in the stress-strain 
curve of a composite lamina: matrix microdamage, matrix macroscopic cracking (modes I 
and II), and axial fiber failure (mode I). Each of these mechanisms can be accommodated 
by partitioning the total dissipated energy density, Ws , into portions associated with each 
mechanism. It is assumed that the evolution of these failure mechanisms (i.e., macroscopic 
cracking and fiber breakage) yields an immediate reduction in the load-carrying capability 
of a local subvolume where the mechanism is active. 

EST in combination with Cohesive Zone methods, such as the discrete cohesive zone 
method (DCZM ),, 20 can be used to account for both in-plane and inter-ply damage. However, 
failure due to delamination is not considered in the second order analysis, since the material 
system studied here did not exhibit significant delamination. 


A. EST Evolution Equations for a Fiber-Reinforced Lamina 

EST is a multiple-internal state variable (ISV) formulation, which utilizes separate ISVs for 
modeling the effects of damage and failure. Four ISVs are used; S, for damage, Sj, for 
fiber failure, S™, for transverse matrix failure, and S'^\ for matrix shear failure. To arrive at 
the evolution equations for the ISVs, the elastic strain energy density must be defined for 
a material, which may contain cohesive cracks. Therefore, the elastic strain energy We is 
comprised of a contribution from the continuum W and any possible cohesive cracks W J M . 
The plane stress, clastic strain energy density in the continuum is defined as, 


kk — + E 2 2 (S)el 2 + Gi 2 ( 5 ) 7 2 2 ) + <512^11622 

where stress in the laminae are related to strain assuming plane stress conditions. 


on — Q 1A11 + Q 12^22 

022 = Q 12^11 + 0,22^22 
T~12 = Q 66712 

where 712 is the engineering shear strain and 
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where Eu is the axial clastic modulus, E 22 is the transverse clastic modulus, z/ 12 is the 
Poisson’s ratio, v 2 \ is the transverse Poisson’s ratio, and G\ 2 is the elastic shear modulus. 
After assuming that the quantity ^ 12^21 << 1, Equations (11) simplify, 


Q 11 = En 
Q 22 = e 22 

Q 12 = ^ 12^22 

Qffi — Gi 2 


( 12 ) 


Note that only the transverse and shear moduli (E 22 and G 12 ) are functions of S since 
matrix microdamage only accrues in the matrix of the laminae. The Poisson’s ratio is 
assumed to evolve such that the quantity Q 12 = E 22 vi 2 remains constant; however, this 
restriction can be relaxed if deemed necessary. The degraded moduli are related to the 
virgin moduli (E 220 and G 120 ) and the ISV through a set of microdamage functions (e s (A) 
and g s (S)). 

E 22 = E 22 oe s (S) (13) 

Gi 2 = Gi 2 og s (S ) (14) 

Degrading E 22 and G i- 2 exclusively is consistent with the intralaminar damage typically 
observed in PMC laminates. 

The elastic strain energy density of the cohesive cracks is defined as the recoverable energy 
per unit crack surface area smeared over the entire element. 
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where, <${ is the fiber direction effective crack opening displacement, 1™ is the matrix trans- 
verse direction effective crack opening displacement, and 5fj, is the matrix shear direction 
effective crack opening displacement. li e+90 ^ is the length of a line running perpendicular 
to fiber direction in the element that intersects two edges of the element and the integration 
point, and le ; is the length of a line that is parallel to the fiber direction in the element 
that intersects two edges of the element and the integration point. The cohesive cracks 
are not active until some initiation criterion is met. Prior to that, the crack tip opening 
displacements S J M are zero, and the do not contribute to the elastic strain energy density. 
The tractions in Equations (15)-(17) can be related to the secant stiffness’ using triangular 
traction-separation laws, k J M shown in figure 7. 
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Hence, the total elastic strain energy density in the continuum is given by 


W-V —v (-E’nCii + E 2 2{S)el 2 + Gi 2 (5)7i 2 2 ) + Qv2^U^22 
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Invoking the stationarity of potential with respect to each ISV, gives the ISV evolution 


equations. 


1 d 63 | 2 ri d'9s "N oq2 

' 220 dS r +ll2Gl20 dS r J- 3Sr 

(22) 

1 dk I zf 2 _ 1 

2 li e+90O) dS{ 1 

(23) 

1 dk? 2 
2 ifdsr ' ~ ' 

(24) 

1 dk]} 2 1 

2 li^dsy} 11 

(25) 


The use of a reduced ISV S r — has been employed in Equation (22). Sicking 21 has shown 
that the use of this reduced ISV yields polynomial forms of the microdamage functions in 
Equations (13) and (14). Thus, Equation (22) becomes a polynomial equation that can be 
readily solved for S r for a given strain state (622,712)- No mixed-mode law is incorporated 
in this work; thus, mode I and mode II failures are uncoupled. Future incrementation 
of EST will implement a mixed-mode fracture law to introduce this coupling. However, 
since the mode I and mode II displacements are uncoupled here, S{,S™, and S 7 /} can be 
determined from Equations (23) - (25), a priori , in terms of the applied strain using the 
traction separation laws and kinematics. 

Using the chain rule and the fact that 
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by traction separation laws, 15 the cohesive secant stiffnesses are determined. The cohesive 
strengths of the material are t{ c (mode I fiber strength), t™ c (mode I matrix strength), and 
fPic ( m ode II matrix strength), and 
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(31) 
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Evaluating the integrals in Equations (29)-(31), while enforcing k J M = 0 when S J M 
results in expressions for k 3 M in terms of 6 3 M . 
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where G IC is the mode I fracture toughness of the fiber, G ™ c is the mode I fracture toughness 
of the matrix, and is the mode II fracture toughness of the matrix. The thermody- 
namically consistent stiffnesses derived in Equations (32)-(34) can also be derived directly 
from the traction-separation laws using geometry. 

Finally, it is assumed that following failure initiation, the strains are related to the crack 
tip opening displacements by 


= i< #+9 "°>£ f, + 6 { 

(35) 

if £22 = if 4 + ST 

(36) 

if 712 = if 7 i C 2 + 2 ST, 
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where ef x , e^, and 7 f 2 are the strains when, the Hashin-Rotem matrix failure criterion, 
involving contributions from both the transverse (622) and shear (712) strains, is satisfied. 


| +(f) 2 = 1 


^22 \ /7i2\ 2 , n 

Yc =1 622 <0 


(38) 


Yt is the transverse lamina failure strain in tension, Yc is the transverse failure lamina strain 
in compression, and Z is the shear failure strain. The fiber failure criterion only involves the 
axial strain en. 

(v;) 2 = 1 e '^° (39) 

where X T is the maximum allowable axial strain of the lamina. The Hashin-Rotem initiation 
criterion is utilized here for its tractability and computational efficiency, but virtually any 
initiation criterion could be used to govern the onset of matrix failure due to macroscopic 
cracking or axial failure due to fiber breakage. 

Equations (35)-(37) imply that the strain in the continuum remains at the values obtained 
when failure initiates, and that any incremental change in the global strain after failure 
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initiation is used wholly to advance the crack tip opening displacement. To account for 
changes in the continuum strain after failure initiates, it can be assumed that the stress 
state in the cracked body is homogenous and the tractions on the crack tip faces are equal 
to the stresses in the continuum. Then, the strains in Equation (21) can be formulated in 
terms of the cohesive secant stiffnesses and the crack tip opening displacement. However, it 
is assumed that the evolution of strain in the continuum is negligible once cohesive cracks 
form. Equations (35)-(37) can be utilized in Equations (32)-(34) to obtain k J M as functions 
of the global strain at an integration point. 
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Once failure initiates, the effects of failure supersede the effects of microdamage and evo- 
lution of S ceases. The cohesive stiffness in a cracked element is calculated using Equations 
(40)-(42) for a given strain state; then, Equations (18)-(20) and (35)-(37) are used to calcu- 
late the tractions on the crack tip faces and the crack tip opening displacement, ft is assumed 
that the stress state in the integration point subvolume of the element is homogenous, and 
the tractions on the crack tip faces are equal to the stresses in the element. Lastly, the axial, 
transverse, and shear moduli of the element can be calculated: 22 
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where E 22 and G* V2 are the degraded transverse and shear moduli, due to microdamage, when 
Equation (38) is satisfied. 
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It should be noted that EST also accounts for kink band formation under axial compres- 
sion. 23-25 As the lamina is loaded, the fibers in the composite rotate by some angle 0, and 
attain the new direction ef, given by the deformation gradient F in the model. 

e f = • e f0 (46) 

where efo is the fiber orientation in the undeformed configuration and A is the stretch ratio. 

A = v/e ro FF T e f0 (47) 

To model the kink band mechanism, all calculations are then executed in the instantaneous 
fiber frame; therefore, fiber rotation can induce large shear strains, 712. As the shear strain 
increases, more damage is induced leading to a reduction in the local tangent shear modu- 
lus. The increase in local shear compliance allows for further increase in the shear strain. 
Under increasing axial compression, this mechanism leads to a runaway (negative feedback) 
instability, and a kink band eventually precipitates. 

In summary, Equations (38) and (39) mark the transition from evolving microdamage 
to failure to macroscopic cracking. Prior to failure initiation, Equation (22) is used to 
calculate the microdamage reduced ISV S r , and the failure ISVs Sj, S™, and Sj} remain 
zero. Equations (13) and (14) are used to calculate the degraded transverse and shear moduli. 
Subsequent to matrix failure initiation, microdamage growth is precluded, and S r remains 
at £*, the value of S r when Equation (38) was satisfied. The degeneration of the transverse 
and shear moduli, resulting from matrix transverse and shear cracking, is calculated using 
Equations (44) and (45). Finally if Equation (39) is satisfied, the axial modulus is calculated 
using Equations (43) as fiber breakage evolves in the element. Once the material moduli 
have been calculated using the appropriate evolution equations, the stresses can be updated 
accordingly using Equations (12). 

V. Material damage and failure characterization 

General procedures for laminates with a cut-outs, like a hole, are shown in £gure(9). 
There are two stages to the analysis, first is the measurement and modeling to obtain input 
parameters from coupon level tests and micromechanics predictions. These are then used 
in the second stage, which is the prediction of OHCS, OHTS and any other progressive 
failure scenario. There are four categories of inputs(£gure 9), geometry, elastic material 
properties, damage properties, and failure properties. Some of the parameters are measured 
while others are derived based on the measured values. Different input details are required 
in each category ranging from constitutive fiber and matrix properties in the micro level to 
the laminate properties at the coupon level as shown in figure 10. 

Due to propriety nature of the composite material used, no information about the material 
will be provided. However, a complete set of EST input and respective test requirements are 
listed in table 1, and are described in subsequent sections. 

A. Measuring lamina in-situ shear response 

In all the methods described , the common and critical data required is the in-situ elastic- 
damage response of the matrix. Characterization of the in-situ lamina shear response is done 
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using the procedure described by Ng et al. 9 Tensile tests of dog-bone shaped specimens of 
±45% were conducted to obtain the in-plane lamina shear response, and the tests followed 
the guidelines of ASTM D3518 and D3039 for strain gage requirements and gage section 
dimensions. Dogbone-shaped specimens were used rather than standard coupon-shaped 
specimens because previous work showed a tendency of specimen failure within the grips for 
straight-sided coupons and within the gage section for dog-bone coupons. 

A speckle pattern was also applied on one side of the specimens to capture full held 
strain data using an optical camera. Subsequently, the speckle images were analyzed using 
the ARAMIS digital image correlation (DIG) software. Strain gage data was compared to, 
and agreed with, the strains obtained through DIG. However, strain gages did not cover 
the full range of strains seen in the experiments and failed prior to sample ultimate failure. 
Hence, DIG measurements were used for calculating full held strain data and used for the 
analysis. 

Using the shear response of the lamina, the in-situ shear response of the damaging matrix, 
which is shown in figure (11a) was extracted as described in. 9 From this, the in-situ matrix 
equivalent stress vs equivalent strain response can be calculated (Ng et al 9 ) and is shown in 
figure (lib). 

B. Damage parameters 

Items 1-3 in table 1 refer to the Schapery damage evolution curves described by the polyno- 
mial forms of e s and g s as 


e s (S r ) — e s o ± e s \S r ± e S 2 ± e s s ± e S iS^ ± e S bS^. (48) 

9s(S r ) = 9s 0 ± 9slS r ± g s 2 Sr ± g s 3&r ± 9s4:S^ ± <? s 5< (49) 

The shear microdamage function g s was obtained from the shear stress-shear strain 
(T 12 - 712 ) response of [45°/-45°] angle-ply coupon tests as recommended by Sicking . 21 The 
transverse tensile micro-damage function e s requires testing of [30°/-30°] angle-ply coupon 
tests, also recommended by Sicking. 21 However, e s can also be obtained from the in-situ ma- 
trix equivalent stress-equivalent strain response. Using the non-linear in-situ matrix proper- 
ties in a cubic close-packed repeating unit cell, and applying transverse loading, the lamina 
transverse response curve can be obtained. The transverse non-linear response curve is then 
used to calculate e s following Sinking’s procedure. 21 The two curves are shown in figure 12a 
& 12b. 

C. Intralaminar failure strain 

Intralaminar failure strains are measured using DIG analysis of tests listed in items 4,6 
and 8 of table 1. Compressive failure strain in the fiber direction can be obtained from 
[0] compression tests, which require thicker coupons such that kink band formation can be 
captured using DIG. Alternatively, compressive failure strain can instead be found using a 
micro- mechanics finite element model (FEM) 1 ' or analytical model 16 of the composite. The 
transverse compression response (item 7) is generally taken to be the same as in tension, 
but the failure strengths are different and the tensile failure strength (or strain) is measured 
from a coupon test. 
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D. Intralaminar failure toughness 

Fiber direction tension failure toughness can be measured from Compact Tension Specimen 
(CTS) 18 or Single Edge Notch Tension (SENT) tests of [90/0] s laminates. The transverse 
direction tension toughness can be obtained from Single Edge Notch Tension (SENT) test. 
However, in this study it is assumed that this value is the same as the interlaminar mode I 
toughness because in both cases, toughness is due to matrix failure only. The transverse di- 
rection compressive toughness is assumed to be the same as the transverse tension toughness. 
Shear direction toughness is difficult to measure because it is difficult to obtain a pure shear 
stable crack propagation. Hence, shear toughness is back calculated from a simulation of 
the SENT test of a [0] ply laminate. Before the shear test and the corresponding simulation 
are conducted, items 1-12 should be completed, so that shear toughness is the last unknown 
variable that needs characterization. 

E. Interlaminar failure toughness 

Interlaminar failure toughness and strength are measured using standard test procedures. 
For mode I toughness and strength, a Double Cantilever Beam (DCB) test is conducted. For 
mode II toughness, an End Notch Flexure (ENF) test is used. Further, for mode II strength, 
a Single Lap Joint (SLJ) test is conducted, with a short intact region so that the join failure 
is strength dominated. 

VI. Finite Element Model and results 

EST was integrated into the Abaqus FEM software using the user material Fortran sub- 
routine UMAT for the implicit solver and VUMAT for the explicit solver. 26 The EST Fortran sub- 
routines are compiled into a static library and linked to Abaqus through the abaqus _v6 . env 
environment hie. At each integration point, during each FEM iteration, the UMAT/VUMAT 
calls the main EST subroutine. 

An open hole laminate is modeled using layered shell elements (S4R) with three integra- 
tion point through the thickness of each lamina. For OHTS simulation, a total displacement 
of 2 mm was applied over lOOOsec using the ^DYNAMIC, IMPLICIT keyword in Abaqus with 
the parameter APPLICATION = QUASI-STATIC. For OHCS, a ramped total displacement of 
2 mm was applied using the ^EXPLICIT, EXPLICIT keyword in Abaqus. For compression 
loading, the explicit solver was used to avoid convergence issues cause by a sudden post- 
peak instability. 

Results from the simulations were compared with experiments for four different layups. 
The nomenclature used to label laminates describe the percentage of 0°,45°&90° plies in 
a laminate. For example laminate (10/70/20) is made of 10% 0° plies, 70% 45° plies and 
20% 90° plies. The maximum load is effect by the stacking sequence. 

Results corresponding to the zeroth order, first order, and second order models are shown 
in figure (13). Although all the predictions are close to experiment, the first and second order 
models are seen to be a better match with experiments (within 5%), especially in laminates 
with high percentage of 0° plies, as shown in figure 13. The results are normalized by the 
compressive strength of the 0° ply with </> = 1° obtained from the Considere construction ie; 
the zero ply un-notched strength. 
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The difference between the zeroth order and second order analysis is more pronounced 
in tension. Results corresponding to the zeroth order and second order models are shown in 
figure (14). The results are normalized by the tensile strength of the 0° ply obtained from an 
experiment ie; the zero ply un-notched tensile strength. Second order analysis results were 
found to be within 5% of the experimental results. First order model for OHTS was not 
pursued, since the tensile failure mechanisms are better captured by second order model. 

VII. Conclusion 

In this paper, different modeling approaches and coupon level test requirements for pre- 
dicting compressive (OHCS) and tensile (OHTS) strength of laminates with a hole are 
outlined. A zeroth order approach which utilized Considere construction and Lekhnitski 
approximations is described. This approach is useful for a efficient approximation of the 
OHCS and OHTS, as it provides conservative values for OHCS and OHTS. 

A first order approach based on global-local micro-mechanical modeling is also outlined. 
This approach can be used for OHCS prediction and is found to be within 5% of experimental 
values. First order approach can capture kink band formation which is the dominant failure 
mechanism in compression. However, this approach does not capture the damage in off-axis 
plies within a laminate. 

A unified tension compression formulation which captures both damage and failure was 
finally described. EST formulation addressed both pre-peak and post-peak non-linearities 
appropriately. Matrix micro-damage, predominantly responsible for the observed pre-peak 
non-linearity in a composite lamina, was accounted for with a single ISV representing the 
dissipated potential associated with micro-damage, along the lines of the original Schapery 
Theory formulation. The relationship between the transverse and shear moduli of the lam- 
ina were related to the ISV through a pair of experimentally-obtainable micro-damage func- 
tions. Three major in-plane failure mechanisms applicable to continuous fiber-reinforced, 
laminated, polymer matrix composites were identified: mode I matrix cracks, mode II ma- 
trix cracks, and fiber breakage. A failure initiation criterion was used to mark the transition 
from a damaging continuum to a damaged continuum with an embedded discontinuity. Cur- 
rently, a quadratic failure criterion is employed, but future work will focus on utilizing a more 
physics-based initiation criterion. After failure initiation, micro-damage evolution ceases and 
separate ISVs are introduced to incorporate the effects of the three major failure mechanisms. 
Evolution of the failure ISVs is based upon traction-separation laws that are material specific 
and that can be measured. Good agreement is seen with first order predictions of OHCS, and 
with second order predictions of both OHCS and OHTS. The material used for this study 
did not show significant delamination. This however is not true for all material systems. 
Future studies will be aimed at incorporating delamination capability using cohesive zone 
elements, as reported in. 20 This ” third order” model will be able to capture all possible 
failure modes associated with PMC laminates, including delamination. 
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Item 

Parameter 

Test 

Damage: Schapery evolution polynomial 

1 

Shear 

[+45 2 / - 45 2 ] s tension 

2 

Trans, tension 

[+30 2 / — 30 2 ] s tension 

3t 

Trans, compression 

[+30 2 / — 30 2 ] s compression 

Intralaminar failure strain 

4 

Fiber tension X et 

[0] tension 

5 

Fiber compression X ec 

[0] compression or FEM 

6 

Trans, tension Y et 

[90] tension 


Trans, compression Y ec 

N/A 

8 

Shear Z e 

[+45 2 / — 45 2 ] tension 

Intralaminar failure toughness 

9 

Fiber tension Gjf_ t 

CTS [90/0 ] s tension 

lot 

Fiber compression G//_ c 

CCS [9 0/0] , s compression 

lit 

Trans, tension Gj t -t 

SENT [90] tension 

12t 

Trans, compression Gu- C 

N/A 

13* 

Shear Gjj s 

SENT [0] tension & FEM 

Interlaminar failure 

14 

Mode I Gj,ai 

DCB 

15 

Mode II G n 

ENF 

16 

Mode II strength r c 

SLJ 

f assumed the same as in tension 


| assumed same as interlaminar failure 


Table 1: Parameters required for modeling 
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kink band 
width 

Figure 1: Typical kink band in a carbon fiber composite. 



Figure 2: Considere construction 
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Figure 3: Infinite plate with a hole under compression 




Lamina: Ortho 
(P-Q) u = 0, v = — S 
(Q-S) u = free, v = free 
(S-R) u = 0, v = 0 
(R-D.A-P) u, v = X symm 


Step 2 


R x 

R„ 



Step 3 


Local liner analysis Local response analysis 


Fiber-ortho, Matrix-iso Fiber-ortho, Matrix-iso 


(A-B) u = A x , v = Ay (A-B) u = A x , t; = A. 


(B-C) u = A x , v = A y 
(C-D) u — A x , v = A y 
( D-A) u = A x , v — A u 


(B-C) F x = R x , F y = R y 
(C-D) u = R x , v = R y 
(D-A) F x = R x , 


Figure 4: Steps in global-local analysis 
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Pinching in micro-region 




Figure 5: Pinching in local analysis when reaction forces are not applied on the free-edge. 
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Figure 6: Typical result in global-local analysis. Top: Load displacement plot from response 
analysis of micro region. Two methods are shown, Standard-dynamic and Riks analysis. 
Middle: Shows the axial stress plots at load points indicated on the load-displacement curve. 
Bottom: shows the kink band formation in the micro region. 
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(b) Mode I matrix fracture. 



(c) Mode II matrix fracture. 

Figure 7: Traction-separation laws used in second order analysis 
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Figure 8: Scaling study result 
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Figure 9: General procedure for OHC & OHT analysis 
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Normalized shear stress 


Geometry 


/ \ 

Micro: fiber dia, mis-alignment. 

Ply: Thickness, vol frac, matrix layer thickness. 

Laminate: Orientations (w.r.t loading), stacking order 


/ \ 

Material 

V y 


Micro: In-situ matrix properties, fiber properties. 

Ply: Linear & Non-linear properties (Hill's anisotropic plasticity) 

Laminate : N/A 

V 


/ \ 

Damage 

v / 


( \ 

Micro: N/A 

Ply: Degradation curves (Schapery curves) 

Laminate : N/A 




Failure 

/ 


( 

Micro: N/A 

Ply: Failure toughness & strength 

Laminate : Inter-ply failure toughness & strength 


Figure 10: Measurements required for a typical laminate analysis 



(a) In-situ shear stress-strain 



Figure 11: Matrix properties 
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Figure 12: Schapery damage curve 



Laminate 

Figure 13: OHCS results 
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Laminate 


Figure 14: OHTS results 
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